Import libraries and data

rm(list=ls())

#define packages desired
dependencies <- c("tidyverse",
                  "plotly",
                  "ggmap",
                  "lubridate",
                  "sp",
                  "rgdal",
                  "geosphere",
                  "rgeos"
                  )

#check if pacakges are installed - load if so, install+load if not)
for (i in dependencies) {
  if (i %in% row.names(installed.packages())){
    eval(bquote(library(.(i))))
    message(paste("loaded package",i))
    }  else {
    install.packages(i)
    eval(bquote(library(.(i))))
  }
}

#read in data
nyc_squirrels <- readr::read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-10-29/nyc_squirrels.csv")

#set mapbox token (free from https://account.mapbox.com/)
#Sys.setenv('MAPBOX_TOKEN'="<SET_YOUR_TOKEN_HERE>")

Data wranglin’

#correct date format
nyc_squirrels <- nyc_squirrels %>% mutate(date = mdy(date))

nyc_squirrels$primary_fur_color <- replace(nyc_squirrels$primary_fur_color, is.na(nyc_squirrels$primary_fur_color), "unknown")

nyc_squirrels$age <- gsub("?", NA, nyc_squirrels$age, fixed=T)

Daily distribution by shift

nyc_squirrels %>% mutate(date=as.factor(date)) %>% group_by(date) %>%
  plot_mapbox(lon=~long, lat=~lat,
              color=~shift,
              text=~paste0("</br>", 
                           location, "</br>",
                           primary_fur_color, "</br>",
                           shift, "</br>"),
              colors = c('gold', 'grey4'),
              mode='markers',
              frame=~date,
              visible=T) %>%

  layout(title="Daily distribution by Shift",
    mapbox=list(
      center=list(lon=~median(long), lat=~median(lat)), 
      zoom =13, 
      bearing=-30,
      pitch=45,
      style='basic'),
    updatemenus = list(
        list(type='buttons',
                direction = "right",
                yanchor = "bottom",
                x = 1,
                y = 0,
                buttons=list(
                  
                  list(method = "relayout",
                      args = list(list(mapbox.style = "basic")),
                      label = "Basic"),
                  
                  list(method = "relayout",
                      args = list(list(mapbox.style = "dark")),
                      label = "Dark"),
                  
                  list(method = "relayout",
                      args = list(list(mapbox.style = "satellite")),
                      label = "Satellite")
                  )
             ))
    ) %>%
  
  animation_opts(
    1000, easing = "elastic", redraw = T
  ) %>%
  
     add_annotations(text = "Hold Ctrl to change pitch and bearing",
                     xref = "paper", yref = "paper",
                     x = 0, xanchor = "left",
                     y = 0, yanchor = "top",
                     legendtitle = TRUE,
                     showarrow = FALSE,
                     font = list(color = '#264E86',
                              family = 'sans serif',
                              size = 14))

Click and drag to pan, hold Ctrl to change camera angle and orientation. This plot shows the distribution of squirrels in the am (gold) and pm (black) surveys on each day. The slider at the bottom can be used to change observational day.

Daily distribution by Primary Fur Color

nyc_squirrels %>% mutate(date=as.factor(date)) %>% group_by(date) %>%
  plot_mapbox(lon=~long, lat=~lat,
              color=~primary_fur_color,
              text=~paste0("</br>", 
                           location, "</br>",
                           primary_fur_color, "</br>",
                           shift, "</br>"),
              colors = c('grey4', 'tan3', 'grey', 'purple'),
              mode='markers',
              frame=~date,
              visible=T) %>%

  layout(title="Daily distribution by Primary Fur Color",
    mapbox=list(
      center=list(lon=~median(long), lat=~median(lat)), 
      zoom =13, 
      bearing=-30,
      pitch=45,
      style='basic'),
    updatemenus = list(
        list(type='buttons',
                direction = "right",
                yanchor = "bottom",
                x = 1,
                y = 0,
                buttons=list(
                  
                  list(method = "relayout",
                      args = list(list(mapbox.style = "basic")),
                      label = "Basic"),
                  
                  list(method = "relayout",
                      args = list(list(mapbox.style = "dark")),
                      label = "Dark"),
                  
                  list(method = "relayout",
                      args = list(list(mapbox.style = "satellite")),
                      label = "Satellite")
                  )
             ))
    ) %>%
  
  animation_opts(
    1000, easing = "elastic", redraw = T
  ) %>%
  
     add_annotations(text = "Hold Ctrl to change pitch and bearing",
                     xref = "paper", yref = "paper",
                     x = 0, xanchor = "left",
                     y = 0, yanchor = "top",
                     legendtitle = TRUE,
                     showarrow = FALSE,
                     font = list(color = '#264E86',
                              family = 'sans serif',
                              size = 14))

Click and drag to pan, hold Ctrl to change camera angle and orientation. This plot shows the distribution of squirrels by primary fur color on each day. The slider at the bottom can be used to change observational day.

Daily distribution by Location

nyc_squirrels %>% mutate(date=as.factor(date)) %>% group_by(date) %>%
  plot_mapbox(lon=~long, lat=~lat,
              color=~location,
              text=~paste0("</br>", 
                           location, "</br>",
                           primary_fur_color, "</br>",
                           shift, "</br>"),
              colors = c('brown3', 'green3', 'grey'),
              mode='markers',
              frame=~date,
              visible=T) %>%

  layout(title="Daily distribution by Location",
    mapbox=list(
      center=list(lon=~median(long), lat=~median(lat)), 
      zoom =13, 
      bearing=-30,
      pitch=45,
      style='basic'),
    updatemenus = list(
        list(type='buttons',
                direction = "right",
                yanchor = "bottom",
                x = 1,
                y = 0,
                buttons=list(
                  
                  list(method = "relayout",
                      args = list(list(mapbox.style = "basic")),
                      label = "Basic"),
                  
                  list(method = "relayout",
                      args = list(list(mapbox.style = "dark")),
                      label = "Dark"),
                  
                  list(method = "relayout",
                      args = list(list(mapbox.style = "satellite")),
                      label = "Satellite")
                  )
             ))
    ) %>%
  
  animation_opts(
    1000, easing = "elastic", redraw = T
  ) %>%
  
     add_annotations(text = "Hold Ctrl to change pitch and bearing",
                     xref = "paper", yref = "paper",
                     x = 0, xanchor = "left",
                     y = 0, yanchor = "top",
                     legendtitle = TRUE,
                     showarrow = FALSE,
                     font = list(color = '#264E86',
                              family = 'sans serif',
                              size = 14))

Click and drag to pan, hold Ctrl to change camera angle and orientation. This plot shows the distribution of squirrels by location on each day. The slider at the bottom can be used to change observational day.

Nearest neighbor data prep

neighboors=NA
for( i in unique(nyc_squirrels$date)){

date <- nyc_squirrels %>% filter(date==i)
  sp.mydata <- date 
  coordinates(sp.mydata) <- ~long+lat
  d <- distm(sp.mydata)
  min.d <- apply(d, 1, function(x) order(x, decreasing=F)[2])
  newdata <- cbind(date, date[min.d,], apply(d, 1, function(x) sort(x, decreasing=F)[2]))
 
  colnames(newdata) <- c(colnames(date), colnames(date), 'distance')
  
  neighboors <- rbind(neighboors,newdata[,c(3,6,8,9,39,44,45,73)])
}

for_neighbor_plot <- as.data.frame(cbind("pair" = c("Same Age", "Both Adults", "Both Juveniles", "Same Fur", "Both Black Fur", "Both Gray Fur", "Both Cinnamon Fur"),
                            "freq" = c(length(which(neighboors$age==neighboors$age.1))/nrow(neighboors),
                            length(which(neighboors$age==neighboors$age.1 & neighboors$age=="Adult"))/length(which(neighboors$age=="Adult")),
                            length(which(neighboors$age==neighboors$age.1 & neighboors$age=="Juvenile"))/length(which(neighboors$age=="Juvenile")),
                            length(which(neighboors$primary_fur_color==neighboors$primary_fur_color.1))/nrow(neighboors),
                            length(which(neighboors$primary_fur_color==neighboors$primary_fur_color.1 & neighboors$primary_fur_color.1=="Black"))/length(which(neighboors$primary_fur_color=="Black")),
                            length(which(neighboors$primary_fur_color==neighboors$primary_fur_color.1 & neighboors$primary_fur_color.1=="Gray"))/length(which(neighboors$primary_fur_color=="Gray")),
                            length(which(neighboors$primary_fur_color==neighboors$primary_fur_color.1 & neighboors$primary_fur_color.1=="Cinnamon"))/length(which(neighboors$primary_fur_color=="Cinnamon"))),
                            
                            "avg.dist"= c(NA, mean(neighboors$distance[which(neighboors$age=="Adult")]),
                                          mean(neighboors$distance[which(neighboors$age=="Juvenile")]),
                                          NA,
                                          mean(neighboors$distance[which(neighboors$primary_fur_color=="Black")]),
                                          mean(neighboors$distance[which(neighboors$primary_fur_color=="Gray")]),
                                          mean(neighboors$distance[which(neighboors$primary_fur_color=="Cinnamon")]))
                                              
                            ))

for_neighbor_plot <- for_neighbor_plot %>% mutate(freq=as.numeric(as.character(freq)))
for_neighbor_plot$pair <-  factor(for_neighbor_plot$pair, levels = c("Same Age", "Both Adults", "Both Juveniles", "Same Fur", "Both Black Fur", "Both Gray Fur", "Both Cinnamon Fur"))

neighboors %>% filter(!is.na(age)) %>%
  plot_ly(y=~distance,
          x=~primary_fur_color,
          type='violin',
          color=~age,
          colors=c('grey4','tan3','grey','purple'),
          text=~paste0("</br>", unique_squirrel_id)
  ) %>%
  layout(
    title='Nearest Neighbor Distance',
    yaxis=list(title='Distance from closest neighbor (m)'),
    xaxis=list(title='Primary Fur Color')
  )

neighboors %>% filter(!is.na(age)) %>%
  plot_ly(y=~distance,
          x=~age,
          type='violin',
          color=~primary_fur_color,
          colors=c('grey4','tan3','grey','purple'),
          text=~paste0("</br>", unique_squirrel_id)
  ) %>%
  layout(
    title='Nearest Neighbor Distance',
    yaxis=list(title='Distance from closest neighbor (m)'),
    xaxis=list(title='Age Class')
  )

Frequency of nearest neighbor sharing characteristics

for_neighbor_plot %>% 
  
  plot_ly() %>%
    add_trace(x=~pair, y=~freq, type = 'bar',
          #marker = list(color=c('green1', 'green2','green3', 'black', 'grey4', 'grey3', 'tan3')),
          text=~round(as.numeric(as.character(avg.dist)),2), textposition='auto') %>%
  
  layout(title="Nearest Neighbor Similarity Frequency",
         yaxis= list(title="Frequency of Occurance", range=c(0,1)),
         xaxis= list(title="")
  
  )

This plot represents the frequency a given squirrels nearest neighbor shared the characteristics listed on the xaxis. Numbers indicate the average nearest neighbor distance in meters for a given category.

nyc_squirrels %>% mutate(date=as.factor(date)) %>% group_by(date, shift) %>% summarise(n=n()) %>%
  plot_ly(labels = ~shift, 
          values = ~n,
          marker=list(colors=c('rgb(255,215,0)','rgb(30,30,30)')),
          frame=~date) %>%
  add_pie(hole = 0.4) %>%
  layout(title = "Daily counts per Shift",  showlegend = T,
         xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
         yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE))

Daily proportion of counts per shift. Gold taken during AM, black during PM

nyc_squirrels %>% mutate(date=as.factor(date)) %>% group_by(date, primary_fur_color) %>% summarise(n=n()) %>%
  plot_ly(labels = ~primary_fur_color, values = ~n, frame=~date, marker=list(colors=c('rgb(30,30,30)','rgb(205,133,63)','rgb(100,100,100)','rgb(160,32,240)'))) %>%
  add_pie(hole = 0.4) %>%
  layout(title = "Daily Proportions of Primary Fur Color",  showlegend = T,
         xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
         yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE))

Daily breakdown of Primary Fur Color